Closed-shell ring coupled cluster doubles theory with range separation applied on weak 

intermolecular interactions 
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We explore different variants of the random phase approximation (RPA) to the correlation energy derived 
from closed-shell ring-diagram approximations to coupled cluster doubles theory. We implement these variants 
in range-separated density-functional theory, i.e. by combining the long-range random phase approximations 
with short-range density-functional approximations. We perform tests on the rare-gas dimers He2, Nei, and 
Ar2, and on the weakly interacting molecular complexes of the S22 set of Jurecka et al. [Phys. Chem. Chem. 
Phys. 8, 1985 (2006)]. The two best variants correspond to the ones originally proposed by Szabo and Ostlund 
[J. Chem. Phys. 67, 4351 (1977)]. With range separation, they reach mean absolute errors on the equilibrium 
interaction energies of the S22 set of about 0.4 kcal/mol, corresponding to mean absolute percentage errors of 
about 4%, with the aug-cc-pVDZ basis set. 
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I. INTRODUCTION 



In the last decade, there has been a revived interest in the 
random phase approximation (RPA) and other related approx- 
imations for calculating the electron correlation energy of 
atomic, molecular and solid-state systems [1^1]. One par- 
ticularly appealing feature of RPA is its correct description of 
dispersion forces at large separation [42-44]. However, RPA 
is a poor approximation to short-range correlations [1], and, 
in a Gaussian localized basis, RPA calculations have a slow 
convergence with respect to the basis size [2]. A promising 
strategy is thus to combine a long-range RPA-type approxima- 
tion with a short-range density-functional approximation [15- 
17, 28, 31, 33], hence avoiding the inaccurate description 
and slow basis-set convergence of short-range correlations in 
RPA. 

Among the different formulations of RPA, the one based 
on a ring-diagram approximation to coupled cluster doubles 
(CCD) theory [14, 45-47] is particularly attractive since it 
avoids the numerical integration over the adiabatic connec- 
tion, and in principle is amenable to a fast algorithm [14]. 
However, due to the fact that the ring approximation breaks 
the antisymmetry property of the coupled-cluster amplitudes, 
several non-equivalent variants of ring CCD can be con- 
structed, especially when the exchange terms are included. In 
this paper, we explore these various ring CCD variants for 
closed-shell systems, and show that some of them correspond 
to the RPA correlation energy expressions originally proposed 
by Szabo and Ostlund [48, 49]. We apply these closed-shell 
ring CCD variants in the context of range-separated density- 
functional theory, and test them on rare-gas dimers and on 
the weakly interacting molecular complexes of the S22 set of 
Jurecka et al. [50]. 



II. THEORY 

We first show how to rigorously combine a long-range 
CCD calculation with a short-range density functional (for 
details on range-separated density-functional theory, see e.g. 
Refs. 33, 51, 52). We start from a self-consistent range- 
separated hybrid (RSH) calculation [52] 



£rsh = mm{(<I)|r + Vext + + £^',,[«o]), 



(1) 



where f is the kinetic energy operator, Vext is the external 
potential operator (e.g., nuclei-electron interaction), Wj^ is a 
long-range electron-electron interaction operator, ^^^^[n] is 
the associated short-range Hartree-exchange-correlation den- 
sity functional, and <I) is a single-determinant wave function 
with density The long-range interaction is constructed 
with the error function, n'J'^(r) - erf(jir)/r, where // is a pa- 
rameter whose inverse gives the range of the separation. The 
minimizing RSH single-determinant wave function is denoted 
by (Do and its associated (approximate) density by hq. In prin- 
ciple, the exact ground-state energy can be obtained from the 
RSH energy by adding the long-range correlation energy E^^ 



E - -Ersh + E 



(2) 



Several formally exact expressions can be derived for E^J. The 
one that is most convenient for applying coupled-cluster the- 
ory is 



+AE- 



](r) An(r)c/r, 



(3) 
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where ^I*'' is the ground-state wave function of the long-range 
interacting HamiltonianH'''[«] - f + W^J^ + Ve^i + V^^^Xn] with 
the short-range Hartree-exchange-correlation potential opera- 
tor ^^^.^[n] = Jv^^.^[n](r)n(r)c/r written with the density op- 
erator n(r) and v^^.^[n](r) = 6E^^^[n]/dn(r). The long-range 
wave function ^P'"^ is associated with the exact density «. In 
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Eq. (3), the last two terms are the variation of the energy func- 
tional, AEfl - [«] - [nn], and the variation of the 

' Hat Hat"- Hxc^ 

associated potential expectation value due to the variation of 
the density from the RSH one to the exact one. An - n — hq. 
The contribution of these last two terms is expected to be small 
since it is of second order in An 

^E'mc- j vS,,[«](r)A«(r)t/r = 

1 rr 6Ef^ [no] 

-o ^ , An{r)An{r')drdr' + 0(An'). (4) 

2 J J dnir)6n{r') 

Using a spin-unrestricted CCD ansatz (see Appendix A for 
a review of CCD theory) for the long-range wave function, 
I^P^cd) = exp(f 2)100), where = {l/4)Zijabitp'alaialaj 
is the cluster operator for double excitations written in terms 
of the long-range amplitudes (ffj')'^ and occupied (/, j) and 
virtual (a, b) RSH spin-orbital creation and annihilation oper- 
ators, we approximate the long-range correlation energy as 

<CCD = mH''[no]\¥l,cu} - mH^noWo). (5) 

In Eq. (5), the variation of the density has been neglected, 
« R! «o (and thus the contribution of Eq. (4) vanishes), which 
seems appropriate if we define the coupled-cluster density as 
the projected one, <cE)o|n(r)|*l'[^cD> = <<I'o|n(r)|<l>o> = «o(r), 
which does not vary since the CCD wave function does not 
contain single excitations. The long-range correlation energy 
can be calculated as, for real spin orbitals. 



<ccD = ^tr[B>'T"] = itr[K'^T'^l, 



(6) 



where BJJ, = ™d K^', ji, - (a^lO)'"^ ^re the matrices 

of antisymmetrized and non-antisymmetrized two-electron in- 
tegrals with long-range interaction n'Jjg(r), respectively, and 
^/■fl jh ~ ^^ij^^'^ ampUtude matrix. The second equality 

in Eq. (6) is due to the antisymmetry property of the coupled- 
cluster amplitudes r'^ with respect to the exchange of the in- 
dices i and / These amplitudes can be determined by the usual 
coupled-cluster equations, replacing the normal Hamiltonian 
by the long-range one W'^lno], which amounts to using the 
RSH orbital eigenvalues and the long-range two-electron inte- 
grals. The present range-separated CCD method can be seen 
a special case of the more general range-separated coupled- 
cluster approach of Goll et al. [53] which also includes single 
excitations and possibly perturbative triples. 

We now consider the ring-diagram approximation for 
closed-shell systems. A number of closed-shell ring CCD 
variants can be defined. In the ring approximation without 
exchange terms, the direct RPA (dRPA, also sometimes re- 
ferred to as RPA or time-dependent Hartree) amplitudes for 
spin-singlet excitations, 'Tjjj^p^, are obtained by the following 
Riccati equation [14] 

Ijy-lr , It Ir Irj-ilr i Ix^i' W Ir.lx^i" ^^^^ \rr^\): f\ /n\ 

+ ^ AdRPA+ *dRPA ^ + *^dRPA ^dRPA-"' y'' 

with the spin-adapted matrices '/t!' ., - 2{ab\ijy'^ and 



Ae, 



iajh 



ia,jb 

where Aeiajb 



(Cfl - ei)6ij6ab 



is the matrix of the RSH orbital eigenvalue differences 
and a, b refer now to occupied and virtual spatial orbitals, re- 
spectively). Contracting the dRPA amplitudes with the non- 
antisymmetrized two-electron integrals 'K'"^ gives the dRPA 
long-range correlation energy (also referred to as dRPA-1 in 
Ref. 54) 



^f.dRPA 



(8) 



Contracting the dRPA amplitudes with the spin-singlet- 
adapted antisymmetrized two-electron integrals 'B'' - 



ia,jb 



2{ab\ijf - {ab\jif gives the dRPA-nSOSEX (or just SOSEX 
for short) long-range correlation energy [23, 28] 



"cSOSEX 



(9) 



Similarly, in the ring approximation with exchange terms, 
that we will refer to as RPAx (also sometimes referred to as 
RPA or time-dependent Hartree-Fock), the singlet and triplet 
amplitudes 'T^ 



'RPAx ^'^RPAx obtained by the equations 



1 1} Ir I 1 A Ir 1 rplr ,1 rplr 1 A Ir i 1 t''" 1 d li" 1 t'1'' — O (^^\\ 
■ ' ^RPAx""" *RPAx ' *RPAx " ^ RPAx " 



and 



3|}lr , 3 i Ir 3rplr 13x11" 3*li-,3iT'lr 3|}lr 3rplr _ ij (^^\ 
D -I- /V ^ RPAx "I" ^RPAx ^ "I" ^RPAx " ^RPAx " \^ ^ ' 

Where 'AJ.;, ., = Aeu,jb+2{ib\ajf' -{ib\jaf\ ^A'^ ., = Ae^a.jb- 
{ib\iay\ and ^B^^,ji, - -{ob\jiy'. Using these amplitudes in 
the CCD correlation energy expression of Eq. (A9) gives what 
we call the RPAx-11 long-range correlation energy (see, also, 
Refs. 48, 49, 54, 56) 

^c,RPAx-II = ['^'' ^T^RPAx + 3^B'' ^T^RPAx] ' (^2) 

which is equivalent to the plasmon formula expression of 
McLachlan and Ball [57]. Using the same amplitudes in the 
alternative CCD correlation energy expression of Eq. (Al 1) 
gives another RPAx correlation energy which is the second 
approximation proposed by Szabo and Ostlund [Eq. (3.22) 
of Ref. 49] as a zeroth iteration of the self-consistent RPA 
scheme [58-60] 



■"c,RPAx-S02 



ltrf'K"'TijpJ. (13) 



Equations (12) and (13) are not equivalent because the ring 
approximation does not preserve the antisymmetry of the am- 
plitudes with respect to the exchange of two spin-orbital in- 
dices. Using the same amplitudes in place of the singlet and 
triplet restricted amplitudes in the CCD correlation energy ex- 
pression of Eq. (A23) gives another alternative RPAx correla- 
tion energy corresponding to the first approximation proposed 
by Szabo and Ostlund [Eq. (3.20) of Ref. 49, or Eq. (17) of 
Ref. 48] which is an alternative zeroth iteration of the self- 
consistent RPA scheme 



■"c,RPAx-S01 



= itr['B"(iT|^p^,-3T|jp^,)]. (14) 
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FIG. 1: (Color online) Interaction energy curves of Hei, Ne2, and ki^ calculated by the full-range (left) and range-separated (right) RPA 
methods with the aug-cc-pV6Z basis set. Cubic splines are used to interpolate between the calculated points. The accurate curves are from 
Ref. 55. 



This last variant is the one preferred by Szabo and Ostlund 
because in a supermolecule approach it consistently gives a 
dispersion coefficient Ce identical to the one given by the 
Casimir-Polder formula applied with the RPAx polarizabil- 
ities of the fragments, which is not the case for the other 
variants RPAx-II and RPAx-S02. On the other hand, among 
the three RPAx methods proposed here, RPAx-S02 has the 
advantage of involving only singlet excitations and thus is 
not subject to triplet instabilities. The RPAx method of 
Refs. 15, 33, that we will rename RPAx-I here, is yet an- 
other alternative correlation energy expression that involves 
only singlet excitations, but for which, as far as we know, the 
numerical integration over the adiabatic connection cannot be 



avoided (although in practice a single-point quadrature works 
well [31]). It can be shown that the SOSEX, RPAx-I, RPAx-II, 
RPAx-SOl and RPAx-S02 correlation energies all correctly 
reduce to the MP2 correlation energy at second order in the 
electron-electron interaction, but dRPA does not. Finally, we 
note that an another RPAx correlation energy variant first pro- 
posed by Fukuda et al. [61] and defined as 2£'c,rpax ii - £c,mp2 
has also been discussed in the literature [41, 47^9]. It ob- 
viously correctly reduces to the MP2 correlation energy at 
second order, but numerical experience [41] shows that this 
variant gives very inaccurate correlation energies. 
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TABLE I: Interaction energies (in kcal/mol) for the complexes of the S22 set from the range- separated RPA methods with the aug-cc-pVDZ 
basis set. For comparison, range-separated CCD results (without the ring approximation) are also reported. The geometries of complexes are 
taken from Ref. 50 and the reference interaction energies in the rightmost column are taken as the CCSD(T)/CBS estimates of Ref. 62. Mean 
errors (ME), mean absolute errors (MAE) and mean absolute percentage errors (MA%E) are given. 



No. Complex dRPA SOSEX RPAx-I RPAx-II RPAx-SOl RPAx-S02 CCD Reference 

Hydrogen-bonded complexes (HB7) 



1 


(NH3)2 


-2.87 


-2.92 


-3.07 


-3.26 


-3.20 


-3.18 


-3.20 


-3.17 


2 


(H20)2 


-5.16 


-5.23 


-5.33 


-5.42 


-5.40 


-5.39 


-5.41 


-5.02 


3 


Formic acid dimer 


-20.30 


-20.55 


-20.81 


-20.98 


-20.94 


-20.98 


-20.94 


-18.80 


4 


Formamide dimer 


-16.51 


-16.68 


-17.03 


-17.48 


-17.32 


-17.27 


-17.25 


-16.12 


5 


Uracil dimer C2;, 


-21.03 


-21.36 


-21.80 


-22.58 


-22.04 


-22.15 


-22.00 


-20.69 


6 


2-pyridoxine/2-aminopyridine 


-17.07 


-17.28 


-17.81 


-18.89 


-18.25 


-18.20 


-18.08 


-17.00 


7 


Adenine/thymine WC 


-16.53 


-16.73 


-17.29 


-18.25 


-17.81 


-17.69 


-17.62 


-16.74 




ME 


-0.28 


-0.46 


-0.80 


-1.33 


-1.06 


-1.05 


-0.99 


0.00 




MAE 


0.42 


0.53 


0.83 


1.33 


1.06 


1.05 


0.99 


0.00 




MA%E 


3.7% 


4.3% 


5.6% 


8.6% 


6.8% 


6.6% 


6.4% 


0.0% 


Complexes with predominant dispersion contribution (WIS) 












8 


(CH4)2 


-0.30 


-0.31 


-0.42 


-0.56 


-0.53 


-0.51 


-0.51 


-0.53 


9 


(C2H4)2 


-0.97 


-1.02 


-1.28 


-1.66 


-1.52 


-1.47 


-1.45 


-1.50 


10 


Benzene/CH4 


-0.92 


-0.98 


-1.23 


-1.75 


-1.47 


-1.43 


-1.40 


-1.45 


11 


Benzene dimer C21, 


-1.27 


-1.38 


-2.05 


-4.28 


-2.72 


-2.61 


-2.40 


-2.62 


12 


Pyrazine dimer 


-2.99 


-3.10 


-3.78 


-6.12 


-4.49 


-4.34 


-4.14 


-4.20 


13 


Uracil dimer C2 


-8.22 


-8.46 


-9.38 


-11.93 


-10.25 


-10.13 


-9.94 


-9.74 


14 


Indole/benzene 


-2.58 


-2.75 


-3.70 


-7.12 


-4.64 


-4.48 


-4.17 


-4.59 


15 


Adenine/thymine stack 


-9.38 


-9.68 


-10.97 


-15.14 


-12.23 


-12.02 


-11.72 


-11.66 




ME 


1.21 


1.08 


0.43 


-1.53 


-0.20 


-0.09 


0.07 


0.00 




MAE 


1.21 


1.08 


0.43 


1.53 


0.20 


0.13 


0.13 


0.00 




MA%E 


34.3% 


31.2% 


13.9% 


31.7% 


3.1% 


2.5% 


3.8% 


0.0% 


Mixed complexes (MI7) 


















16 


Ethene/ethyne 


-1.31 


-1.36 


-1.48 


-1.67 


-1.58 


-1.57 


-1.55 


-1.51 


17 


Benzene/H20 


-2.90 


-2.96 


-3.16 


-3.52 


-3.34 


-3.30 


-3.29 


-3.29 


18 


Benzene/NH3 


-1.83 


-1.88 


-2.11 


-2.57 


-2.33 


-2.29 


-2.27 


-2.32 


19 


Benzene/HCN 


-4.20 


-4.31 


-4.54 


-4.98 


-4.72 


-4.71 


-4.65 


-4.55 


20 


Benzene dimer C21. 


-1.92 


-2.00 


-2.39 


-3.40 


-2.77 


-2.70 


-2.61 


-2.71 


21 


Indole/benzene T-shape 


-4.54 


-4.65 


-5.17 


-6.57 


-5.66 


-5.57 


-5.44 


-5.62 


22 


Phenol dimer 


-6.48 


-6.62 


-7.07 


-8.16 


-7.49 


-7.43 


-7.35 


-7.09 




ME 


0.56 


0.47 


0.17 


-0.54 


-0.11 


-0.07 


-0.01 


0.00 




MAE 


0.56 


0.47 


0.17 


0.54 


0.11 


0.09 


0.11 


0.00 




MA%E 


15.8% 


13.5% 


5.0% 


13.6% 


2.7% 


2.2% 


2.6% 


0.0% 




total ME 


0.53 


0.40 


-0.04 


-1.15 


-0.44 


-0.39 


-0.30 


0.00 




total MAE 


0.75 


0.71 


0.47 


1.15 


0.44 


0.41 


0.40 


0.00 




total MA%E 


18.7% 


17.0% 


8.4% 


18.6% 


4.1% 


3.7% 


4.3% 


0.0% 



III. COMPUTATIONAL DETAILS 

All calculations have been done with a development version 
of MOLPRO 2008 [63], implementing equations (7)-(14). We 
first perform a self-consistent RSH calculation with the short- 
range Perdew-Burke-Ernzerhof (PBE) exchange-correlation 
functional of Ref. 64 and add the long-range RPA correlation 
energies calculated with RSH orbitals. The range separation 
parameter is taken at jj - 0.5 bohr"', according to previous 
studies [65], without trying to readjust it. For the rare-gas 
dimers, we also carry out full-range RPA calculations using 
PBE orbitals [66] for comparison. The Riccati equation (7) 
is solved by decomposing the matrix 'L'' into diagonal and 
off-diagonal parts and iteratively extracting 'Tjjj^p^ from its 



product with the diagonal part and updating it in the other 
terms, and similarly for Eqs. (10) and (11). For RPAx-I calcu- 
lations, the adiabatic-connection integration is performed by 
a 8-point Gauss-Legendre quadrature for the rare-gas dimers, 
and by a single-point quadrature [Eq. (14) of Ref. 31] for the 
S22 set. We use the correlation-consistent basis sets of Dun- 
ning [67, 68]. Core electrons are kept frozen (i.e. only exci- 
tations of valence electrons are considered). Basis set super- 
position error (BSSE) is removed by the counterpoise method. 
The geometries of the complexes of the S22 set are taken from 
Ref. 50. The geometries of the isolated monomers are fixed 
to those in the dimers, thus the so-called monomer deforma- 
tion energy is not included in the interaction energy. For each 
method, mean error (ME), mean absolute error (MAE) and 
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mean absolute percentage error (MA%E) are given using as 
a reference the CCSD(T) values extrapolated to the complete 
basis set (CBS) limit of Takatani et al. [62]. 

In our present, most basic implementation, the computa- 
tional cost of all the RPA methods used here formally scales as 
N^NI for large basis sets, where A^,, and A^,, are the numbers of 
virtual and occupied orbitals, respectively. The computational 
cost of the CCD (or CCSD) method without the ring approxi- 
mation is higher and it scales as N^.N^ for large basis sets [69]. 
Of course, far better scalings should be obtained by using 
integral-direct methods and resolution-of-identity/Cholesky- 
decomposition techniques [14]. 



IV. RESULTS 

The interaction energy curves of Hea, Nei, and Ar2 cal- 
culated by the full-range and range-separated RPA methods 
are compared in Fig. 1. We use the large aug-cc-pV6Z ba- 
sis set to ensure that the full-range calculations are converged. 
Full-range dRPA and SOSEX strongly underestimate the in- 
teraction energies, while full-range RPAx-II and RPAx-SOl 
strongly overestimate them. The best full-range methods are 
RPAx-I and RPAx-S02, which is in agreement with the re- 
cent study of HeBelmann [41, 70]. In passing, we note that 
the full-range RPAx-I method better performs for Ne2 and 
Ar2 when using PBE orbitals than when using HE orbitals, as 
done in Ref. 33. Range separation greatly improves the accu- 
racy of all the RPA variants. However, range-separated dRPA 
and SOSEX still underestimate the interaction energies, and 
range-separated RPAx-II significantly overestimates the inter- 
action energy of Ar2. Range-separated RPAx-I, RPAx-SOl, 
and RPAx-S02 give the most reasonable interaction energy 
curves. 

The interaction energies for the complexes of the S22 set 
calculated with the range-separated RPA methods with the 
aug-cc-pVDZ basis set are given in Table I. For comparison, 
range-separated CCD results (without the ring approximation) 
are also reported. Although the aug-cc-pVDZ basis set may 
appear small, range-separated RPA methods are weakly de- 
pendent on the basis size [15, 33], and indeed it was estimated 
in Ref. 3 1 that when going from the aug-cc-p VDZ to the aug- 
cc-pVTZ basis set the range-separated RPAx-I interaction en- 
ergies of the S22 set are lower by at most 7%, and the corre- 
sponding total MA%E decreases by less than 2%. Therefore, 
we believe that the aug-cc-p VDZ basis set is sufficient to com- 
pare the different range-separated RPA methods. 

The S22 set includes seven hydrogen-bonded complexes 
(HB7 subset), eight weakly interacting complexes with pre- 
dominant dispersion contributions (WIS subset), and seven 
mixed complexes featuring also multipole interactions (MI7 
subset). The trends are quite different for the HB7 subset on 
the one hand, and the WIS and MI7 subsets on the other hand. 
It was previously argued that the general overestimation of the 
interaction energies of hydrogen-bonded complexes is due to 
the approximate short-range density functional [31, 71]. The 
fact that dRPA and SOSEX give the smallest MAEs for the 
HB7 subset is thus not believed to be significant but rather due 



to a compensation of errors between an underestimated long- 
range contribution and an overestimated short-range contri- 
bution. This is corroborated by the relatively large overes- 
timation of the interaction energies of this subset by range- 
separated CCD which should most accurately describe the 
long-range correlation energies. We will thus focus our anal- 
ysis on the WIS and MI7 subsets. 

For the WIS and MI7 subsets, dRPA gives largely under- 
estimated interaction energies, with MA%Es of 34.3% and 
15.S%, respectively. SOSEX barely improves dRPA with 
MA%Es of 31.2% and 13.5%, respectively. This may not be 
surprising since, in the limit of large separation, SOSEX only 
adds exponentially decaying exchange interactions between 
the monomers, but does not change the coupled-cluster am- 
plitudes and thus does not change the polarizabilities of the 
monomers. The RPAx-I method of Refs. 15, 33, which in- 
corporates exchange effects in the monomers, greatly reduces 
the underestimation of the interaction energies, with ]VlA%Es 
of 13.9% and 5.0%, respectively. The RPAx-II variant, which 
may be seen as the most straightforward way of defining a 
closed-shell ring CCD with exchange terms, is disappoint- 
ingly inaccurate. It overestimates the interaction energies by 
about the same amount that dRPA underestimates them. Fi- 
nally, the two variants RPAx-SOl and RPAx-S02 give re- 
markably accurate interaction energies, with MA%Es of 3. 1 % 
and 2.7% for RPAx-SOI, and 2.5% and 2.2% for RPAx-S02. 
They are globally as accurate as range-separated CCD without 
the ring approximation. However, it must be noted that RPAx- 
SOl and RPAx-S02 tend to overestimate dispersion energies, 
while RPAx-I underestimates them. Therefore, increasing the 
basis size will likely increase the MA%Es of RPAx-SOl and 
RPAX-S02, while it will decrease the MA%E of RPAx-I. 



V. CONCLUSION 

We have studied various RPA variants that can be cast in 
the form of closed-shell ring CCD approximations. We have 
tested these variants with range separation, i.e. by com- 
bining a long-range RPA-type approximation with a short- 
range density-functional approximation, on rare-gas dimers 
and on the weakly interacting complexes of the S22 set. 
Among all these variants, the ones first proposed by Szabo 
and Ostlund [4S, 49], called here RPAx-SOl [Eq. (14)] and 
RPAx-S02 [Eq. (13)], give the most accurate dispersion ener- 
gies. The other variants tend to either strongly underestimate 
(dRPA and SOSEX) or strongly overestimate (RPAx-II) the 
interaction energies. For comparison, we have also reported 
results from the RPAx-I method of Refs. 15, 33, which is not 
based on a ring CCD approximation but on the adiabatic con- 
nection formula, and which gives reasonable interaction en- 
ergies as well. From a practical point of view, RPAx-S02 
appears to be the most convenient variant since, contrary to 
RPAx-I, it does not use any numerical adiabatic-connection 
integration and, contrary to RPAx-SOl, it involves only sin- 
glet excitations and is thus not subject to triplet instabilities. 
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Appendix A: CCD correlation energy 

In this appendix, we review several equivalent CCD corre- 
lation energy expressions, in view of justifying the different 
ring CCD variants. 



1. CCD correlation energy in spin-orbital basis 

The spin-unrestricted CCD wave function ansatz is 

l*PccD> = exp(f2)|cD>, (Al) 

where lO) is a single-determinant reference wave function, 
and T2 is the cluster operator for double excitations which is 
written in a spin-orbital basis as 



^2 = ^2 fi'lKciihlaj, 



(A2) 



ijab 



where /, j and a, b refer to occupied and virtual spin-orbitals. 



respectively, and the amplitudes f.j must be antisymmetric 
with respect to any exchange of two indices: f^'^ - -fjf - 



fba 



fha 



-t"j' — t'J. . The CCD correlation energy is obtained by the 
transition formula 



-CCD 



- ^2<«^l|/;< = ^tr[BT], 



(A3) 



ijab 



where Biajb - {ab\\ij} are the antisymmetrized two-electron 
integrals over real spin orbitals and Tiajb - t"j is the ampli- 
tude matrix. Using the antisymmetry of the amplitudes, the 
CCD correlation energy can also be written as 

= i2<fl/7|/;)f;;^ = ltr[KT], (A4) 

ijab 

where Kjajb - {ab\ij} are the two-electron integrals. 



C = 



TT.TT 


Ctt.u 





^ 


U.TT 


iUi 























CiT.n 


CiT.iT . 



(A5) 



and can be brought to a block-diagonal spin-adapted matrix 
C = C U by the orthogonal transformation 



1 1 ^ 
1-10 
11 
1-1 



(A6) 



Applying this transformation to the matrix B gives the follow- 
ing decomposition into singlet and triplet excitations 



B 



/''BOO ^ 
^B 
^B 

-^B 



(A7) 



where Biajb - '2{ab\ij) - {ab\ji) and Biaj, - -{ab\ji}, with 
j referring now to occupied spatial orbitals and a, b to virtual 
spatial orbitals. Notice the minus sign for the last triplet block. 
Using Kramers symmetry for spin-conserving real coupled- 
cluster amplitudes (see, e.g., Ref. 72), one can show that spin 
adaptation of the matrix T leads to a similar form 

'TOO ^ 






-^T j 



(A8) 



where ^Tiaj, - Ti1a1J■lb^ + TiiaTjibi and ^Tiajb = T^/toTjT'-'T ~ 
TiMjm = Ti^aijib-f. The CCD correlation energy can thus be 
expressed as 



itrf'B'T + S^B^Tl. 



(A9) 



Spin adaptation of the matrix K gives only a contribution from 
the singlet excitations 



K = 



r 'K 0^ 









(AlO) 



where ^Kiajb - 2{ab\ij}, which leads to an alternative form 
for the CCD correlation energy 



(All) 



2. CCD correlation energy in spatial-orbital basis for 
closed-shell systems 

a. Expression in terms of the singlet and triplet amplitudes 

For spin-restricted closed-shell calculations, all the matri- 
ces in the spin-orbital excitation basis encountered so far (e.g.. 



b. Expression in terms of the restricted amplitudes 

In practice, the CCD correlation energy is normally calcu- 
lated starting from the spin-restricted closed-shell CCD wave 
function ansatz 

l*l'ccD) = exp('^f2)|cD>, (A12) 
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where the restricted cluster operator is written in a spatial- 
orbital basis as 



ijab 



(A13) 



constructed with the singlet single-excitation operator 



'«f = 4(«It«'T+«IA) = ^4. (A19) 



where = fl^^fl,| -H al^an is the singlet excitation opera- and faihj is the triplet double-excitation operator 



tor and t^j are the restricted amplitudes which must be sym- 
metric with respect to the exchange of both i,j and a,b, i.e., 
^fl^ = ^t'^", but not antisymmetric with respect to the ex- 
change of only two indices. The CCD correlation energy is 
obtained by the transition formula 



(A20) 



constructed with the triplet single-excitation operators 



7CCD 



= (OIZ/I^Fccd) - {mm = (OIH^falO) 

- 2 (2<«^l'7) - {ab\ji)) X = tr ['B '^t] (A14) 



ijab 



where i^riaji = 



ab 



(A21a) 



(A21b) 



c. Expression in terms of the singlet and triplet restricted 
amplitudes 

Another equivalent correlation energy expression can be 
obtained by decomposing the restricted amplitudes into spin- 
singlet and spin-triplet components. Indeed, the restricted 
cluster operator can be decomposed as (see, e.g., Ref. 73, 74) 

ijab 



where ^-^f'j are singlet restricted amplitudes 



U ji U ' 



(A16) 



which are totally symmetric (i.e., ^'^t"j - ^'^fjf - '-^f*" - 
and ^.Rfa* ^re the triplet restricted amplitudes 



3,R^fl6 _ R.ab _ R.ab 
U Ji ij ' 



(A17) 



which are totally antisymmetric (i.e., • t"j - - ■ t^f - 
_XRfbj ^ XRfip, In Eq. (A15), S aibj is the singlet double- 
excitation operator 



'Jaibj - J „i J i,j - -t^ail^bj. 



(A18) 



f : — a ,an. 



(A21c) 



Using the symmetry properties of '•'^rj* and "^ '^f"*, it is easy 
to check that Eqs. (A13) and (A15) are equivalent. Combin- 
ing Eqs. (A16) and (A17) leads to the decomposition of the 
restricted amplitudes into spin components 

^^j-^i^^tfj-'^^f), (A22) 



and the CCD correlation energy [Eq. (A14)] can thus be writ- 
ten as 



(A23) 



This corresponds to the definition of singlet and triplet con- 
tributions to the correlation energy, - (l/2)tr['B '•'^T] and 
EJ - -(l/2)tr['B ''''*T]. By using the symmetry properties 
of '•'^T and ■'•'^T, one can show that they are equivalent to 
the more usual expressions in terms of the restricted ampli- 
tudes (see, e.g., Ref. 75): E^^ = (l/4)tr[(iB - 3^B)'*T] and 
El = (3/4)tr[(iBH- WT]. 
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